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28 Abstract 

29 A land surface model’s ability to simulate states (e.g., soil moisture) and fluxes (e.g., 

30 runoff) is limited by uncertainties in meteorological forcing and parameter inputs as well as 

3 1 inadequacies in model physics. In this study, anomalies of terrestrial water storage (TWS) 

32 observed by the Gravity Recovery and Climate Experiment (GRACE) satellite mission were 

33 assimilated into the NASA Catchment land surface model in western and central Europe for a 7- 

34 year period, using a previously developed ensemble Kalman smoother. GRACE data 

35 assimilation led to improved runoff correlations with gauge data in 17 out of 18 hydrological 

36 basins, even in basins smaller than the effective resolution of GRACE. Improvements in root 

37 zone soil moisture were less conclusive, partly due to the shortness of the in situ data record. In 

38 addition to improving temporal correlations, GRACE data assimilation also reduced increasing 

39 trends in simulated monthly TWS and runoff associated with increasing rates of precipitation. 

40 GRACE assimilated root zone soil moisture and TWS fields exhibited significant changes in 

4 1 their dryness rankings relative to those without data assimilation, suggesting that GRACE data 

42 assimilation could have a substantial impact on drought monitoring. Signals of drought in 

43 GRACE TWS correlated well with MODIS Nonnalized Difference Vegetation Index (NDVI) 

44 data in most areas. Although they detected the same droughts during warm seasons, drought 

45 signatures in GRACE derived TWS exhibited greater persistence than those in NDVI throughout 

46 all seasons, in part due to limitations associated with the seasonality of vegetation. 


47 
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48 1 . Introduction 

49 Seasonal and interannual variability in terrestrial water storage (TWS) is of critical 

50 interest in water resource analysis and seasonal hydrological forecasts because TWS — which 

5 1 includes soil moisture, groundwater, surface water and snow — is an important hydrological 

52 indicator in its own right: volume of water stored in snowpack or groundwater, for example, 

53 reflects present hydrological conditions and can be used to infer the potential for future 

54 hydrological stress. TWS is also important because of its role in other aspects of the 

55 hydrological cycle. Its status can affect infiltration rates and subsurface flow, with associated 

56 impacts on runoff and recharge rates. TWS anomalies can also affect the hydrological cycle 

57 through soil moisture feedbacks on the atmosphere. One of the important aspects of TWS is its 

58 unique dynamics. Soil moisture and groundwater are low-pass filters on the terrestrial 

59 hydrological cycle that gradually remove high frequency variability associated with atmospheric 

60 forcing as depth increases (Eltahir and Yeh, 1999; Wu et al. 2002). This dynamic means that 

61 TWS acts as a “memory” component of the terrestrial hydrological cycle, with implications for 

62 land-atmosphere interactions (Koster and Suarez, 2001) and predictability in certain regions 

63 (Dinneyer 2000; Dirmeyer et al., 2009; Koster et al., 2000b; Koster et al., 2010a). 

64 Interactions among components of TWS not only re-distribute water spatially but also 

65 increase the complexity of the hydrological cycle. Groundwater, which accounts for a major part 

66 of TWS (Rodell and Famiglietti, 2001; Rodell et al., 2007; Yeh et al. 2006), can contribute 

67 substantially to stream flow in wet climates (Eltahir and Yeh, 1999). This connection, combined 

68 with the long memory of groundwater variability, means that accurate information on 

69 groundwater can contribute significant skills to seasonal river discharge forecasts (Birkens and 

70 Van Beek, 2009). Groundwater can also move upward to increase soil wetness through capillary 
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71 lift or act as a sink to receive excess soil moisture from the land surface (Schaller and Fan, 2009). 

72 As appreciation for these processes has grown, an increasing number of land surface models 

73 have been developed to account for the impact of groundwater on near surface processes (e.g., 

74 Koster et ah, 2000a; Niu et ah, 2007; Miguez-Macho et ah, 2007; Yeh and Eltahir, 2005). 

75 Including groundwater in a land surface model enables a more complete simulation of the 

76 terrestrial water cycle, but it also subjects the modeled states to additional uncertainties 

77 associated with the added physical processes and parameters. For instance, due to lack of global- 

78 scale groundwater measurements, most models depend on calibration to obtain the temporal 

79 variability and dynamic range of groundwater tables, which may not represent the interactions 

80 realistically, especially under extreme wet or dry conditions. 

8 1 Precipitation data sets are a major source of uncertainty for land surface modeling, and 

82 their impacts on modeled states and fluxes may differ depending on seasons and climates (Fekete 

83 et al., 2004; Gottschalck et al., 2005). Great uncertainty also exists in model physics such as 

84 surface runoff algorithms which are often derived from empirical relationships (Koster et al., 

85 2000a; Niu et al., 2005; Schaake et al., 1996). Stream flow is governed in varying degrees by 

86 topography, rainfall intensity, and soil wetness, making it a difficult process to simulate 

87 efficiently. Due to differences in model physics and parameter values, estimates by various land 

88 surface models exhibit large discrepancies even when models are run using identical forcing data 

89 (Mitchell et al., 2004). The combination of uncertainties in forcing, input parameters and model 

90 physics has led to dramatically different predictions for runoff trends in response to future 

91 climate changes (Hoerling et al., 2009). 

92 The ambiguity in model estimates also complicates drought monitoring, which 

93 increasingly relies on model estimated soil moisture due to the current lack of accurate global 
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soil moisture measurements (Mo, 2008). Although Koster et al. (2010b) provided a more 
optimistic assessment on soil moisture estimates by various models, Mo (2008) indicated that 
while drought indices derived from different models show stronger correlation in the eastern US, 
their correlation is so low in the western US that model based drought indices cannot be used for 
drought monitoring. Drought monitoring is also complicated by the interaction between soil 
moisture and groundwater. Through numerical simulations, Peters et al. (2005) showed that 
groundwater can provide moisture to reduce the impact of short-term droughts, but due to its 
long recovery time groundwater will also act to lengthen and increase the frequency of droughts. 
The importance of groundwater for drought monitoring has been recognized (Houborg et al., 
2011; Svoboda et al., 2002) and efforts are underway to combine information about groundwater 
variability as well as surface vegetation conditions with model estimated soil moisture to form 
comprehensive drought indices (http://www.drought.unl.edu/dm/monitor.html) . Nevertheless, 
such efforts are hindered by the lack of systematic groundwater measurements at continental 
scales, in addition to lack of accurate model based soil moisture estimates. 

In order to capture the unique characteristics of TWS and reduce the uncertainty in model 
estimates, observations are needed to nudge model output towards reality. The GRACE satellite 
system detects temporal water storage changes in the entire vertical profile, including snow 
mass, surface water, vegetation, soil moisture and groundwater (Tapley et al., 2004). It is the 
only remote sensing platfonn that provides consistent monitoring of the Earth’s terrestrial water 
storage, including groundwater. Recognizing the potential for GRACE data to improve the 
simulation of land surface processes, Zaitchik et al. (2008) developed an ensemble Kalman 
smoother (EnKS) to assimilate GRACE into the NASA Catchment model in the Mississippi 
basin, with promising results. The EnKS provides a systematic and dynamic way to disaggregate 
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117 GRACE-derived TWS anomaly estimates into snow, soil moisture, and groundwater 

118 components, so that the simulation of each component of TWS can be positively influenced. 

119 In this study, the EnKS and the Catchment model are applied in western and central 

120 Europe where climate and hydrological conditions differ significantly from the Mississippi area 

121 studied by Zaitchik et al. (2008). As droughts are common in Europe, the unique ability of 

122 GRACE TWS to detect droughts and its potential for drought monitoring are considered in some 

123 detail. The paper is organized as follows: Sections 2 and 3 describe the study domain, ground 

124 based validation data and the land surface model. Section 4 briefly outlines the EnKS method 

125 and filter parameters. Section 5 presents the model simulation results and comparisons with 

126 independent datasets. Comparisons of anomalies of GRACE TWS with those of MODIS NDVI 

127 are also presented. Section 6 concludes with a summary and discussion. 

128 2. Experiment site, GRACE and validation data 

129 Figure 1 shows the simulation domain in western and central Europe. For GRACE data 

130 assimilation, major hydrological watersheds were combined into nine major “basins” at the scale 

131 of GRACE observations, to accommodate the spatial resolution of GRACE TWS, which is about 

132 150,000 km 2 at best (Rowlands et al., 2005; Swenson et al., 2006). Table 1 lists the area of these 

133 basins, ranging from 300,000 to 800,000 km". Several islands and peninsulas such as Great 

134 Britain and Sweden/Norway were not included because GRACE TWS yielded much smaller 

135 dynamic ranges than model estimates, possibly due to the interference of ocean signals. 

136 GRACE TWS used in this study were processed by University of Texas Center for Space 

137 Research (CSR, Release CSR RL04) using a Gaussian filter with a 300 km smoothing radius to 

138 remove the stripes seen in the spherical harmonic coefficient fields (Swenson and Wahr, 2006). 
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139 The anomalies of GRACE TWS were obtained by removing the temporal mean of the gravity 

140 field (including the solid earth and the atmosphere) in 2003-2007 and converted to equivalent 

141 water heights. The 1° gridded GRACE TWS anomalies were mapped to the nine major basins 

142 using area- weighted averaging, and these values were converted to absolute TWS by adding the 

143 2003 - 2007 mean TWS from an open loop (no data assimilation) integration of the model. 

144 Figure 1 also shows the locations of in situ measurements used for validating data 

145 assimilation results, including 18 stream flow stations along three major rivers (Danube, Elbe 

146 and Rhine) and 12 soil moisture sites from the Soil Moisture Observing System - Meteorological 

147 Automatic Network Integrated Application (SMOSMANIA, Calvet et ah, 2007) project. The 

148 streamflow stations (station ids and drainage areas are given in Table 2) were chosen from 

149 Global Runoff Data Center (GRDC) for their length of records. Soil moisture measurements 

150 (started in 2007) are taken at 5, 10, 20 and 30 cm depths and every 30 minutes using impedance 

151 probes. Monthly averaged stream flow and root zone soil moisture (vertically integrated using 

152 the four layer measurements) were used to validate model simulation results. 

153 3. The Catchment model and forcing data 

154 The NASA Catchment model was developed for global scale coupled land/atmosphere 

155 modeling (Koster et ah, 2000a). It simulates water and energy balances on catchment tiles, with 

156 some catchments split by a 1.0°xl.25° atmospheric grid. For the study domain, which consists of 

157 nearly 6000 tiles, the average tile size is around 1500 km". To increase sub-grid heterogeneity, 

158 each catchment contains dynamically changing saturated, transpiring and wilting areas where 

159 different runoff and ET schemes are applied. The model contains three subsurface states for 

160 water balance calculation: surface excess (sfEx) and root zone excess (rtzEx), representing the 
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161 excessive soil moisture relative to the hydrostatic state for the top 2 cm and 100 cm of soils, 

162 respectively, and catchment deficit (catDef) defined as the amount of water (kg/m', averaged 

163 over the catchment) needed to bring the catchment to saturation (assuming sfEx and rtzEx are 

164 zero). Although groundwater is not explicitly simulated, its behavior, i.e., its two dimensional 

165 distribution and associated flow rates, is directly diagnosed from the catDef variable. The model 

166 also has three snow layers for modeling snow water equivalent (SWE) and snow depth. Thus, 

167 modeled TWS can be determined from sfEx, rtzEx, catDef and SWE in conjunction with model 

168 parameters. Lakes and reservoirs are not directly included in simulated TWS because, over large 

169 scales at mid-latitudes, they only constitute a very small fraction of observed TWS variability 

170 (Rodell and Famiglietti, 2001). The impact of GRACE data assimilation on runoff is exerted 

171 through its relationship with modeled states: sfEx, rtzEx, catDef and SWE. 

172 Forcing fields were provided by the Global Land Data Assimilation System (GLDAS, 

173 Rodell et al. 2004). They are based on meteorological fields (temperature, humidity, wind speed 

174 and pressure) obtained from the NASA Global Modeling and Assimilation Office GEOS data 

175 assimilation system (Bloom et al., 2005), radiation fields from the U.S. Air Force Weather 

176 Agency, and precipitation prepared by spatially and temporally downscaling the 2.5° x2. 5°, 5 -day 

177 NOAA Climate Prediction Center (CPC) Merged Analysis of Precipitation (CMAP; Xie and 

178 Arkin, 1997). This GLDAS forcing data set, which has been used in previous data assimilation 

179 experiments (Reichle et al., 2007; Zaitchik et al., 2008), has a 3 hour temporal interval and a 2° x 

180 2.5° spatial resolution. 

181 A few adjustment and corrections were made in this study regarding the Catchment 

182 model and forcing fields. Zaitchik et al. (2008) found that Catchment sometimes does not 

183 provide a large enough dynamic range to match that of GRACE TWS. The same situation was 
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184 observed in this study region as well. To mitigate this deficiency, following Houborg et al. 

185 (2011), the bedrock depth used for the model was uniformly increased by 2 m, which increased 

186 the dynamic range of catDef. To partially compensate for the increase in bedrock depth, a lower 

187 value of the decay factor for saturated conductivity was used for the base flow calculation 

188 (Duchame et ah, 2000). Longwave and shortwave radiation fields were further bias corrected 

189 based on NASA/GEWEX Surface Radiation Budget (SRB, Release-3.0) data by matching their 

190 spatial (for entire simulation area) and temporal averaged means with those of SRB. The goal of 

191 these adjustments and corrections was to achieve reasonable estimates of fluxes (ET and runoff). 

192 Simulations were carried out from August 2002 to July 2009, which is the available 

193 GRACE data period at the start of this study. Since previous forcing data were not available, the 

194 model was first run through 2002 to 2009 and then spun up for 10 years using the forcing fields 

195 from 2002. A different initial condition, based on averaged model states from 2002-2009 on 

196 January 1 which yielded wetter soil moisture conditions than the one mentioned above, was also 

197 tested and the results (including runoff and soil moisture evaluations) were very similar to those 

198 presented here. 

199 4. GRACE data assimilation method 

200 Zaitchik et al. (2008) presented a detailed description of the ensemble Kalman Smoother 

201 (EnKS) developed specifically for assimilating GRACE TWS into the Catchment model. A brief 

202 outline of this assimilation method is presented here. Like an ensemble Kalman filter (EnKF), 

203 the EnKS consists of two steps: forecast and update. In the forecast step, the ensemble of the 

204 model runs forward in time with perturbations added to the states and forcing fields: 
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X‘ T - =M(X i (Ti)+ ,F i ,G) (1) 

205 where M is the model; F represents all the forcing fields and G represents all the static 

206 parameters; T is the time; superscripts (-) and (+) refer to results for the forecast and update, 

207 respectively; X is the vector containing updated states (rtzEx, catDef and SWE) for each 

208 catchment tile, and the superscript i indicates the i th member of the ensemble. srfEx was not 

209 updated in the EnKS because of its very weak correlation with monthly TWS but was included in 

210 model simulated TWS for accuracy. Based on equation (1), the ensemble update equation can be 

211 written as: 

X T : = X T J + K t (Y t - H(X t -‘ )) (2) 

212 where K is the ensemble gain matrix; Y represents observations (GRACE TWS) and H is the 

213 observation operator that converts predicted states to the observation. 

214 The underscores in equation (2) indicate monthly TWS (observed or simulated) averaged 

215 for each major basin because the EnKS used here assimilates temporally integrated observations. 

216 To accommodate the monthly averaged nature of GRACE observations, the EnKS collects 

217 Catchment model predictions of TWS on a first pass through each simulated month (three 

218 collections per month, to mimic GRACE overpass characteristics), calculates the update at the 

219 end of the month, and then iterates through the month a second time, uniformly (for each state) 

220 applying increments to each daily value of model states for each ensemble member. Thus, X 

22 1 (without the underscore) in equation (2) represents daily estimates of model states on each 

222 catchment tile in month T. 
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223 All perturbation parameters and schemes were the same as Zaitchik et al. (2008) and 

224 Reichle et al. (2007), except that an observation (GRACE) error of 15 mm was used here, which 

225 is the average of the two GRACE errors (10 mm and 20 mm) tested by Zaitchik et al. (2008). 

226 5. Results 

227 Two model integrations were performed from August 2002 to July 2009: the open loop 

228 (OL) representing the model-only performance and the data assimilation (DA) with GRACE data 

229 assimilation using the EnKS outlined above. Since GRACE derived TWS values are anomalies 

230 only, simulation results were evaluated using time series correlations with in situ measurements. 

231 5.1 TWS 

232 Figure 2 presents the time series of daily simulated TWS and GRACE monthly 

233 observations for the nine major basins. The open loop run generally captured the seasonal 

234 variability and dynamic range of GRACE TWS. OL differs from GRACE mostly in interannual 

235 variability, especially in Finland, Loire/Seine and Rhone/Po where OL exhibits a marked 

236 increase in TWS in the later modeling period. While data assimilation checked that increase 

237 effectively, consistent with GRACE observations, it failed to reduce simulated TWS to the levels 

238 observed by GRACE in the Upper Danube in 2007 and 2008. This failure was possibly caused 

239 by negligible ensemble spread during dry conditions due to a lack of precipitation to perturb and 

240 the fact that direct perturbations to sfEx and catDef are small. Increasing the direct perturbations 

241 may enable TWS to go lower, but it may also lead to ensemble bias. Nevertheless, in most cases 

242 EnKS was effective in nudging the simulated TWS toward GRACE TWS. 

243 To investigate the cause of the significant increase in OL TWS seen in the Finland, 


244 Vistula, Loire/Seine, and Rhone/Po basins, which was not observed as dramatically in the 
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245 GRACE observations (Figure 2), GLDAS/CMAP precipitation was compared with 1° x 1° 

246 Global Precipitation Climatology Project (GPCP) precipitation data (Adler et ah, 2003) mapped 

247 to the major basins following the same approach as that for GRACE. Figure 3 shows the 

248 comparison of annual (from August to July) precipitation totals in each basin. In general, 

249 GLDAS/CMAP has a negative bias against GPCP in all basins except Turkey. CMAP's low bias 

250 relative to other precipitation products stems from the fact that it does not correct for gauge 

25 1 under-catch (e.g., Yin et ah, 2004). More importantly, the annual variations of GPCP and 

252 GLDAS/CMAP precipitation are well correlated, and both products indicate that precipitation in 

253 the four basins named above increased towards the end of the simulation period. Given that 

254 GRACE TWS also increased in those basins but to a lesser extent, we infer that either: (i) the 

255 model should have retained less water in the land and increased evapotranspiration (ET) and/or 

256 runoff instead, or (ii) the precipitation and GRACE datasets are inconsistent, due to errors in one 

257 or both. 

258 The rate of long-term TWS changes can be more clearly illustrated using the slope of 

259 monthly TWS calculated using Sen’s method (Helsel and Hirsch, 1992; Sen, 1968) as shown in 

260 Figure 4. Slopes with a 0. 1 significance level were identified using the Mann-Kendal test 

261 (Helsel and Hirsch, 1992) and marked in bold symbols. These two methods have been widely 

262 used in analyzing trends in hydro-meteorological data sets (Mishra and Cherkauer, 2010; 

263 Lettenmaier et ah, 1994; Yue and Wang, 2002). Figure 4 shows that the slope of TWS (modeled 

264 or observed) generally becomes smaller as the basin moves from north to south, which resembles 

265 the increasing rate of annual precipitation in each basin (Figure 3), suggesting the strong 

266 correlation of TWS with long-term precipitation. OL TWS generally exhibits larger rates of 

267 increase than GRACE-derived TWS, especially in Finland, Vistula, Loire/Seine and Upper 
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268 Danube, where larger increasing rates of precipitation were observed in the later modeling period 

269 (Figure 3). 

270 5.2 Stream flow and soil moisture 

27 1 Since stream flow is a product of upland surface runoff and subsurface runoff over a 

272 large area, gauged stream flow data are often used to evaluate model performances (Mishra and 

273 Cherkauer, 2010). For the same reason, stream flow measurements were used here to not only 

274 evaluate the impact of GRACE data assimilation on runoff but also provide overall assessment of 

275 the EnKS. Figure 5 shows the correlation of monthly simulated runoff with GRDC gauged data. 

276 Since Catchment does not have a routing scheme, the simulated stream flow is simply a spatial- 

277 aggregation of tile-space (individual land element) runoff over the drainage area. This is 

278 justifiable for monthly statistics, especially in smaller basins where the runoff response time is 

279 less than a month. GRACE data assimilation improved the correlation in all watersheds but one 

280 (D5), with more improvements observed in larger basins along Danube. Improvements in 

281 watersheds such as R6-R1 1, El and E2 (Table 2) with drainage areas smaller than their 

282 corresponding major basins (the scale at which GRACE TWS was generated) indicate that 

283 assimilation of GRACE TWS can influence simulation of land surface processes at sub- 

284 observation scale. The improvements shown in Figure 5 by DA all exceeded the 0.05 

285 significance level based on the William-Hotelling t-test (Steiger, 1980; Van Sickle, 2005). It 

286 should be pointed out that many of the stream flow observations are not independent because 

287 they were measured at various points along the same river. 

288 Improvements in runoff correlations are attributed to the close relationship between base 

289 flow and catDef, which is the model state most affected by assimilation of GRACE TWS. To 
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290 illustrate this, Figure 6 shows the time series of simulated runoff in comparison with GRDC 

291 measurements in Lower Danube 6742800, a sub-basin of the Lower Danube major basin. DA 

292 significantly increased the runoff in the earlier period in accordance with changes in TWS, which 

293 helped improving the overall correlation and also lowered the increasing trend of runoff. Figure 

294 7 shows the trend of runoff by OL, DA and GRDC gauge data in all GRDC basins. Similar to 

295 TWS, model estimates (OL) show higher trends than observed runoff with significant trends 

296 detected for most basins while observed runoff shows no significant trend in any basin. DA 

297 reduced trends in all basins, but did not change the significance of most trends. 

298 An important role of the EnKS is to disaggregate GRACE so that each TWS component 

299 can be nudged towards its true state. To evaluate the vertical disaggregation, correlations of 

300 monthly root zone soil moisture estimated by OL and DA were calculated against in situ 

301 measurements from the SMOSMANIA sites and are given in Table 3. The statistics were 

302 calculated using in situ point data and model estimates at the tile containing the station. GRACE 

303 data assimilation generally did not have a significant impact on monthly correlations of soil 

304 moisture as the correlation of DA is not significantly different from OL at the 0.10 significance 

305 level, except at site URG. The coarser spatial representation of the model and the GRACE data 

306 may not capture the localized nature of station measurements. To alleviate the horizontal scale 

307 mismatch and obtain an overall impact on the entire SMOSMANIA area (about 4000 km“), the 

308 area averaged statistics for OL and DA were also calculated against the averaged in situ 

309 measurements and are given in Table 3 (last row) which shows that GRACE data assimilation 

310 did not change the correlation of averaged soil moisture time series in the sampling area. The 

311 shorter SMOSMANIA data period (3 1 months) makes these statistics less conclusive because the 

312 confidence intervals are very large. 
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313 

3 14 5.3 Water budget 

315 As hypothesized in section 5.1, elevated TWS by OL in Finland and Loire/Seine in the 

316 later modeling period were likely caused by either an underestimation of runoff and/or ET when 

317 precipitation rates increased or by improper increase in the precipitation rates themselves 

318 assuming GRACE data are accurate. When GRACE data assimilation reduced TWS in these 

319 basins, it also decreased ET and runoff estimates because of their positive correlations with 

320 TWS. As a result, the water budget of OL was not preserved by the simulation with GRACE 

32 1 data assimilation. Figure 8 features the annual (August to July) mass imbalance, defined as 

322 simulated water budget (sum of total fluxes and net change in TWS) minus precipitation, of OL 

323 and DA. As expected, OL has nearly zero mass imbalances throughout the entire period and in 

324 all basins while GRACE data assimilation disrupted the water budget, more so in Finland, 

325 Vistula, Loire/Seine and Rhone/Po, despite improving the simulation of TWS (assuming 

326 GRACE data are accurate). Since GLDAS precipitation generally has a negative bias against 

327 GPCP (Figure 3), positive imbalances (i.e., larger ET and runoff) would be preferable to the 

328 negative ones produced by GRACE data assimilation in this case. Unintended impacts of data 

329 assimilation on the water budget are always a danger, demanding the development of creative 

330 new assimilation techniques (e.g., Li et ah, 2011; Pan and Wood, 2006; Zaitchik and Rodell, 

331 2009). 

332 5.4 Drought analysis 

333 Droughts are common in Europe, and several episodes of severe droughts, including the 

334 2003 drought (associated with the 2003 European heat wave, Rebetez et ah, 2006; Zaitchik et ah, 
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335 2006) that spread across western and central Europe and the 2007/2008 droughts that affected 

336 southern and southwestern Europe (SOER Synthesis, 2010), were detected by GRACE TWS 

337 (Figure 2). Because droughts can be defined in a variety of ways depending on what indicators 

338 are taken into account, it can be instructive to compare a new drought observation with a more 

339 common indicator. Here we compare GRACE based TWS with monthly Normalized Difference 

340 Vegetation Index (NDVI) as recorded by the Moderate Resolution Imaging Spectroradiometer 

341 (MODIS) instrument on NASA’s Terra satellite. NDVI is strongly correlated with green 

342 biomass (Tucker, 1979), and is often used in satellite based drought-monitoring (e.g., Brown et 

343 ah, 2008). Basin averaged NDVI was derived by averaging the Level-3 0.05° MODIS NDVI 

344 monthly product (lpdaac.usgs.gov) across the same basins that were used to extract GRACE 

345 observations. 

346 Figure 9 shows the averaged dryness ranks of NDVI and GRACE TWS in the summer 

347 season (April to September) for the 2003 to 2008 period (2002 and 2009 were excluded due to 

348 their incomplete summer seasons). To give equal weight to all monthly rankings, the averaged 

349 ranks in Figure 9 were obtained by first ranking each data set for each month and then averaging 

350 the ranks of summer months. GRACE TWS indicated 2003 as the driest condition in all basins 

35 1 except Loire/Seine, Lower Danube and Turkey, while NDVI only shows 2003 as the most severe 

352 drought in Rhone/Po, Upper Danube and Dnieper and a drought condition in Rhine/Elbe/Oder 

353 and Loire/Seine. The 2007/2008 droughts along the south and southwestern region (in 

354 Rhone/Po, Lower/upper Danube, Dnieper and Turkey) were indicated by both types of 

355 observations. The largest discrepancies between the two sources are in Finland and Vistula 

356 where, despite the increasing trend in precipitation (Figure 3), NDVI shows decreasing trends. 
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357 This is probably due to the fact that vegetation growth in the high latitude and high elevation 

358 regions is limited by energy availability, not by water availability (Karnieli et ah, 2010). 

359 Note that GRACE TWS characterized the 2003 drought in Loire/Seine as less severe than 

360 the 2005 drought (SOER Synthesis, 2010). According to GRACE, the land was very wet in 

361 early 2003 (Figure 2), and as a result dry meteorological conditions took longer to severely 

362 impact total TWS. In this situation, the effect of drought is less evident in the TWS anomaly 

363 than it is in the maximum decline of GRACE TWS from its early spring peak to the lowest value 

364 in the fall, which roughly measures the amount of water lost in the wann season. As seen in 

365 Figure 9, Loire/Seine and Upper Danube, which were at or near the center of the heat wave, 

366 experienced the most significant water loss in 2003. This is one of the advantages using a 

367 physical-based variable for drought monitoring because drought conditions can be evaluated 

368 from other aspects than anomalies. 

369 The reason that we only compared the dryness rank of GRACE and NDVI during the 

370 warm season in Figure 9 is that NDVI is insensitive to water shortage when vegetation is 

371 senescent or when coverage is low (Karnieli et ah, 2010). This can be seen in Figure 10 where 

372 the seasonal cycles of GRACE TWS and NDVI in the Lower Danube basin are presented. 

373 GRACE TWS shows signs of stress in 2007 very consistently over all seasons, in contrast with 

374 NDVI which indicated vegetation stress only after June. GRACE-derived TWS also exhibits 

375 large inter-annual variability and larger dynamic ranges that can provide more infonnation on 

376 drought severity. These qualities, true in most areas (Rodell, 2011), are important both for 

377 drought monitoring and for early detection of drought onset and therefore make GRACE a useful 

378 complement to high-resolution NDVI-based measures of drought, especially in regions with low 

379 vegetation cover or where water is not a limiting factor for vegetation growth. 
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380 Figures 9 and 10 show the dryness ranks based on GRACE TWS data alone. To 

381 demonstrate the potential value of integrating GRACE and other data with a land surface model 

382 for drought monitoring, Figure 1 1 plots the dryness ranks (among 2002 to 2009) of OL and DA 

383 estimated root zone soil moisture (upper panels), which is of particular interest for monitoring 

384 agricultural droughts, and TWS (lower panels), which is an indicator of water depletion in the 

385 deeper subsurface, for November 2007. GRACE DA intensified the drought condition in 

386 Loire/Seine and Upper Danube relative to the open loop. The updates in both the root zone soil 

387 moisture and TWS demonstrate that data assimilation makes it possible to apply GRACE for 

388 monitoring both agricultural and hydrological droughts, and to do so with much greater spatial 

389 resolutions than with GRACE alone. 

390 6. Summary and Discussions 

391 This study demonstrated the value of GRACE TWS for correcting errors in model 

392 estimated TWS and its influence on related land surface processes. In particular, assimilation 

393 significantly improved runoff correlation in most basins, which attests to the overall robustness 

394 of the assimilation technique and the usefulness of GRACE TWS for runoff estimation. The 

395 improved runoff correlation in small watersheds also shows the potential of GRACE TWS to 

396 contribute to simulation of finer scale hydrological processes through data assimilation based 

397 downscaling. Assimilation of GRACE TWS did not improve the correlation of soil moisture 

398 with in situ measurements, perhaps due to the short in situ data record or insufficient spatial 

399 sampling. Although groundwater was not validated directly due to lack of in situ measurements, 

400 the improvements in stream flow estimates suggest more realistic estimates of subsurface water 

401 storage which controls baseflow. 
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GRACE data assimilation had a significant impact on reducing trends of modeled TWS 
and runoff. The original inconsistency between the GRACE and OL trends is caused by 
deficiencies in either the model’s physics, the forcing data or the GRACE data themselves. The 
case presented here represents a relatively short period during which annual precipitation 
increased at a much higher rate in several basins than long tenn annual precipitation trends 
(Mishra and Cherkauer, 2010; Solomon et ah, 2007). The fact that GRACE TWS was able to 
change the trend in runoff suggests that GRACE TWS data, if independently validated, may 
assist in model and forcing evaluation and calibration, which is an important part of climate 
prediction (Mishra and Cherkauer, 2010), especially in regions with scarce observation data. 
However, only models able to simulate groundwater storage can take full advantage of GRACE, 
because assimilation of GRACE TWS requires an analogous model state. 

Monitoring of droughts has suffered from lack of reliable information on the water stored 
below the uppermost soil layer. Since GRACE measures the water storage changes in the entire 
profile, it provides valuable infonnation on drought development beyond what can be seen at the 
surface. Its large dynamic range and inter-annual variability also provides better quantification 
of the severity of water depletion in the subsurface. The continued monitoring of dry conditions 
throughout all seasons, which cannot be achieved using vegetation based indicators, may further 
assist in tracking prolonged droughts and/or providing early signs of drought development. 

While data assimilation of GRACE TWS helps to fill the need for regional to global scale 
infonnation on deep moisture storage variability, it also presents some challenges. Since drought 
indices are derived based on the long tenn climatology of a given variable (Mo, 2008) and the 
GRACE observation period is not long enough to generate its own climatology, GRACE based 
drought indices must be linked to a model simulation that begins well prior to data assimilation. 
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425 This requires that the estimates from GRACE assimilation have the same dynamic range as 

426 GRACE, so that the anomalies from the assimilation period are comparable to the climatology. 

427 To accomplish this, it may be necessary to modify parameters such as the bedrock depth, which 

428 controls the amount of water available from storage to be lost during a drought (Houborg et ah, 

429 2011). The changing trends in DA TWS, as found in this study, may also reduce the dynamic 

430 range and the magnitude of anomalies and thus present a new challenge. Statistical techniques 

43 1 such as cumulative distribution function matching may also be used to ensure that the modeled 

432 and observed climatologies are consistent prior to generating drought indices (Houborg et ah, 

433 2011). Nevertheless, these challenges should not discourage the use of GRACE data 

434 assimilation for drought monitoring because the dryness infonnation provided by GRACE TWS 

435 can lead to more objective and reliable drought indices (Rodell, 2011). 

436 Water budget imbalance caused by GRACE data assimilation is an important issue for 

437 future research because existing flux biases may be exacerbated (assuming precipitation forcing 

438 data were accurately estimated). In this example, we speculate, without the benefit of ET and 

439 runoff observations in Finland and Loire/Seine regions, that a low bias in modeled ET and runoff 

440 might have caused the TWS anomaly to be elevated, which, when corrected by GRACE data 

441 assimilation, further reduced ET and/or runoff. This water budget imbalance might have been 

442 avoided, if observations of ET and runoff were available and assimilated simultaneously with 

443 GRACE TWS. Given that ET and runoff observations are rarely assimilated into land surface 

444 models, a more likely solution would be to remove excess TWS during the assimilation process 

445 in conjunction with increasing simulated ET and/or runoff. Exploring creative new data 

446 assimilation strategies such as this is recommended so that the benefits of GRACE DA can be 
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447 realized while avoiding detrimental effects on modeled water budgets (Li et ah, 2011; Pan and 

448 Wood, 2006; Zaitchik and Rodell, 2009). 

449 
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625 Figures 

626 Figure 1: Study area and major basin boundaries. The blue cross and red triangle represent 

627 locations of GRDC stream flow and SMOSMANIA soil moisture sites, respectively. Numbers 1 

628 to 9 represent the nine major basins given in Table 1 . 

629 Figure 2: Time series of daily simulated TWS and monthly GRACE TWS in the nine major 

630 basins. 

63 1 Figure 3: Comparisons of annual GLDAS and GPCP precipitation in the nine basins. 

632 Figure 4: Slopes of trend for monthly TWS in the nine major basins. Trends with a 0. 1 

633 significance level are marked with bold symbols. 

634 Figure 5: Correlations of monthly simulated runoff with GRDC stream flow. All improvements 

635 by DA exceed the 0.05 significance level. Station ids are given in Table 2. 

636 Figure 6: Monthly time series of estimated runoff in comparison with GRDC gauge data. 

637 Figure 7: Slopes of trend for monthly runoff at GRDC stations. Trends with a 0. 1 significance 

638 level are marked with bold symbols. 

639 Figure 8: Annual mass imbalance (simulated water budget minus precipitation) for OL and DA 

640 in the nine major basins. 

641 Figure 9: Averaged dryness ranks of NDVI and GRACE TWS for the summer growing season 

642 (April to September) during the 2003 to 2008 period and maximum GRACE TWS declines from 

643 spring to fall in each year. 

644 Figure 10: Seasonal cycles of GRACE TWS and NDVI in Lower Danube. 
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645 Figure 1 1 : Dryness ranks of simulated root zone soil moisture and TWS for November 2007 in 

646 the 2002 to 2009 period. 

647 
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650 Figure 1. Study area and major basin boundaries. The blue cross and red triangle represent 

65 1 locations of GRDC stream flow and SMOSMANIA soil moisture sites, respectively. Numbers 1 

652 to 9 represent the nine major basins given in Table 1 . 


653 


33 


654 

655 

656 



date 



date 




date 



date 


date 



1000 


Rhine/Elbe/Oder 


1600 - 
1400 ■ 


co 

S 1200 • 


1000 ■ 


'/yvvV^ 


— i 1 1 1 1 r 

2003 2005 2007 2009 

date 



date 


1000 



date 


Figure 2. Time series of daily simulated TWS and monthly GRACE TWS in the nine major 
basins. 
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659 Figure 3. Comparisons of annual GLDAS and GPCP precipitation in the nine basins. 
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661 Figure 4. Slopes of trend for monthly TWS time series in the nine major basins. Trends with a 

662 0.1 significance level are marked with bold symbols. 
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665 Figure 5. Correlations of monthly simulated runoff with GRDC stream flow. All improvements 

666 by DA exceed the 0.05 significance level. Station ids are given in Table 2. 
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669 Figure 6. Monthly time series of estimated runoff in comparison with GRDC gauge data. 
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674 Figure 7. Slopes of trend for monthly runoff at GRDC stations. Trends with a 0.1 significance 

675 level are marked with bold symbols. 
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678 Figure 8. Annual mass imbalance (simulated water budget minus precipitation) for OL and DA 


679 in the nine major basins. 












40 


680 


681 


Finland 



100 | 
£ 


200 J3 
o 

V 

T) 

■ 300 OT 


year 

Loire/Seine 



• 100 


■ 200 


u 

C 
•«— < 

o 

0) 

T3 


300 V} 
H 


year 
Upper Danube 



100 6 
£ 
<D 

200 5 

73 

0) 

73 

■ 300 CO 


year 


Vistula 



100 § 
£ 
v 

■ 200 J3 

73 
0) 
T* 

■ 300 CO 


year 


Rhone/Po 



100 


• 200 


300 


<D 

5 

o 

<D 

73 

CO 

► 

H 


year 


Dnieper 


"3 

c0 

u 



a 

^a 

0> 

a 

0 

<D 

73 

1 


year 


Rhine/Elbe/Oder 



year 
Lower Danube 



year 


Turkey 



year 


682 Figure 9. Averaged dryness ranks of NDVI and GRACE TWS for the summer growing season 

683 (April to September) during the 2003 to 2008 period and maximum GRACE TWS declines from 

684 spring to fall in each year. 
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688 


689 Figure 10. Seasonal cycles of GRACE TWS and NDVI in Lower Danube. 
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691 Figure 11. Dryness ranks of simulated root zone soil moisture and TWS for November 2007. 
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693 Table 1 : Major basins and their drainage areas. 

694 Table 2: GRDC stations, drainage areas and record lengths. 

695 Table 3: Correlations of monthly averaged simulated soil moisture with observations at 

696 SMOSMANIA sites. Except for the URG site, the OL and DA correlation values are not 

697 significantly different at the 0.10 significance level. 

698 
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699 Table 1 : Major basins and their drainage areas. 

700 


701 

702 

703 

704 

705 

706 

707 

708 

709 

710 


basin ID 

basin name 

Area 

(km 2 ) 

1 

Finland 

498000 

2 

Vistula 

547000 

3 

Rhine/Elbe/Oder 

797000 

4 

Loire/Seine 

393000 

5 

Rhone/Po 

319000 

6 

Lower Danube 

503000 

7 

Upper Danube 

490000 

8 

Dnieper 

721000 

9 

Turkey 

403000 


711 
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712 Table 2: GRDC stations, drainage areas and record lengths. 


713 


714 

Station ID 

GRDC 

number 

drainage area 
(km 2 ) 

record length 
(months) 



Rhine 

715 

R1 

6435060 

160800 

65 


R2 

6335020 

159300 

77 

716 

R3 

6335050 

147600 

77 


R4 

6335060 

144200 

77 

717 

R5 

6335070 

139500 

77 

R6 

6335100 

103500 

77 

718 

R7 

6335150 

98200 

77 

R8 

6335180 

68800 

77 


R9 

6335170 

53100 

77 

719 

RIO 

6335200 

50200 

77 


Rll 

6335400 

34600 

77 

720 


Elbe 

El 

6340110 

131900 

77 

721 

E2 

6340150 

123500 

77 


Danube 


D1 

6742900 

807000 

77 

722 

D2 

6742800 

709100 

84 


D3 

6742500 

658400 

84 

723 

D4 

6742201 

570900 

77 


D5 

6242501 

101500 

53 


724 


725 

726 


727 
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728 Table 3: Correlations of monthly simulated soil moisture with observations at SMOSMANIA 

729 sites. Except for the URG site, the OL and DA correlation values are not significantly different 

730 at the 0.10 significance level. 
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